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ABSTRACT 


Context. High-frequency very-long-baseline interferometry (VLBI) observations can now resolve the event-horizon-scale emission 
from sources in the immediate vicinity of nearby supermassive black holes. Future space-VLBI observations will access highly lensed 
features of black hole images—photon rings—that will provide particularly sharp probes of strong-field gravity. 

Aims. Focusing on the particular case of the supermassive black hole M87*, our goal is to explore a wide variety of accretion flows 
onto a Kerr black hole and to understand their corresponding images and visibilities. We are particularly interested in the visibility on 
baselines to space, which encodes the photon ring shape and whose measurement could provide a stringent test of the Kerr hypothesis. 
Methods. We develop a fully analytical model of stationary, axisymmetric accretion flows with a variable disk thickness and a matter 
four-velocity that can smoothly interpolate between purely azimuthal rotation and purely radial infall. To determine the observational 
appearance of such flows, we numerically integrate the general-relativistic radiative transfer equation in the Kerr spacetime, taking 
care to include the effects of thermal synchrotron emission and absorption. We then Fourier transform the resulting images and analyze 
their visibility amplitudes along the directions parallel and orthogonal to the black hole spin projected on the observer sky. 

Results. Our images generically display a “wedding cake” structure composed of discrete, narrow photon rings (n = 1, 2,...) stacked 
on top of broader primary emission that surrounds a central brightness depression of model-dependent size. At 230 GHz, the n = 1 
ring is always visible, but the n = 2 ring is sometimes suppressed due to absorption. At 345 GHz, the medium is optically thinner and 
the n = 2 ring displays clear signatures in both the image and visibility domains. We also examine the thermal synchrotron emissivity 
in the equatorial plane and show that it exhibits an exponential dependence on radius for the preferred M87* parameters. 
Conclusions. The “black hole shadow” is a model-dependent phenomenon—even for diffuse, optically thin sources—and should not 
be regarded as a generic prediction of general relativity. Observations at 345 GHz are promising for future space- VLBI measurements 
of the photon ring shape, since at this frequency the signal of the n = 2 ring persists despite the disk thickness and non-zero absorption 
featured in our models. Future work is needed to investigate whether this conclusion holds in a larger variety of reasonable models. 


Key words. Physical data and processes: Gravitation — Accretion, accretion discs — Black hole physics — Relativistic processes — 


©ESO 2022 


Galaxies: individual: M87 


1. Introduction 


In 2019, the Event Horizon Telescope collaboration released 
1.3 mm interferometric observations of the supermassive black 
hole M87* at the center of the galaxy Messier 87 (EHT; EHT 
L1), achieving an effective angular resolution comparable to the 
black hole size. These observations revealed the presence of a 
bright ring of approximately 40 yas in diameter that surrounds a 
much darker central region. While these basic image features are 
undisputed and have in fact been independently confirmed (Ar- 
ras et al. 2022; Carilli & Thyagarajan 2022; Lockhart & Gralla 
2022), several aspects of their theoretical interpretation remain 
open. For example, should the dark area be associated with the 
“black hole shadow” (Falcke et al. 2000), as originally proposed 
(EHT L1), or with the apparent position of the equatorial event 
horizon, as the data now seem to suggest (Chael et al. 2021)? 


A plausible range of observational appearances for M87* is 
depicted in Fig. 1. Under the currently favored assumption of 
optically thin emission concentrated very near the event horizon 
(e.g., EHT L5; EHT L8), the characteristic appearance ranges 


between two extremes: a narrow photon ring surrounding a dark 
region inside the critical curve [a “shadow”—see Falcke et al. 
(2000)], and a series of discrete photons rings stacked on top 
of broader emission outside the apparent equatorial horizon [a 
“wedding cake’”—see Gralla et al. (2019, Fig. 1) and Johnson 
et al. (2020, Fig. 3)]. Models of spherically symmetric, infalling 
matter lead to shadows, while models with equatorial, orbiting 
matter generate wedding cakes. General-relativistic magneto- 
hydrodynamic (GRMHD) models generally favor the wedding 
cake over the shadow (Johnson et al. 2020; Chael et al. 2021), 
but the debate has yet to be settled (Bronzwaer & Falcke 2021). 


The purpose of this paper is to more fully flesh out the model 
space between the shadow and wedding cake extremes. Having 
a broad range of models is crucial not just for understanding the 
source but also for accurately forecasting its potential for future 
observations. We are especially motivated by the promise of pho- 
ton ring (orbiting light) measurements on long interferometric 
baselines (Johnson et al. 2020; Gralla 2020), which could pro- 
vide a stringent test of the Kerr hypothesis (Gralla et al. 2020; 
Wielgus 2021). A recent proposal to measure the precise shape 
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of the n = 2 photon ring (formed by light that executes a full 
orbit around the black hole before escaping to a detector) relied 
exclusively on a purely equatorial and perfectly absorption-free 
model of the emission (Gralla et al. (2020)—henceforth GLM 
20). It is imperative to check whether this exciting prospect sur- 
vives the introduction of geometrical thickness, absorption, and 
other potential complications of the real astrophysical flow (see 
also Paugnat et al. 2022). 


In this paper, we construct a set of semi-analytic models that 
bridge the gap between the shadow and wedding cake extremes, 
while also including a “realistic” level of absorption. For select 
density, temperature, and velocity profiles, we compute the syn- 
chrotron emission and absorption based on the assumption of 
a uniformly magnetized disk. Tuning parameters to match the 
total 230 GHz horizon-scale flux density reported by EHT, we 
perform radiative transport to determine the observational ap- 
pearance at both 230 GHz and 345 GHz, the planned frequency 
of future observations (Doeleman et al. 2019; EHT L2). Finally, 
we Fourier transform the images and check whether the photon 
ring signatures are observable on moderate and long baselines. 
We consider both infalling and orbiting matter, as well as geo- 
metrically thin, thick, and fully spherical emission regions. We 
fix the observer inclination to a value of 160°, deemed to be very 
likely for M87* (Walker et al. 2018), while varying the black 
hole spin from near-zero to near-maximal spin. 


Our results have both theoretical and practical implications. 
On the theory side, we confirm that the size of the central bright- 
ness depression is highly model-dependent, while the presence 
of discrete photon rings is generic (Gralla et al. 2019; Johnson 
et al. 2020). We also confirm that the classic “black hole shadow” 
(a sharp intensity drop inside the critical curve) is caused by ex- 
treme special-relativistic redshifting (Narayan et al. 2019; Gralla 
2021), depending only indirectly on general-relativistic effects. 
Furthermore, we find that the intensity drop occurs only in fully 
spherical, infalling models, whereas the dark area is highly dis- 
tinct from the critical curve even for very thick disks assumed to 
contain purely infalling matter. That is, far from being a generic 
prediction for optically thin flows (as claimed, e.g., in Psaltis 
2019; Narayan et al. 2019), the appearance of a dark “shadow” 
filling the interior of the critical curve is in fact a remarkably fine- 
tuned phenomenon. This provides further evidence that the EHT 
observations should not be regarded as the black hole shadow in 
the strict sense given by Falcke et al. (2000). 


On the practical side, our results indicate that the n = 2 pho- 
ton ring is only marginally observable at 230 GHz. For strongly 
magnetized disks, we reproduce previous results (Johnson et al. 
2020; Chael et al. 2021) showing a prominent n = 2 photon ring, 
but we find that in the weak-magnetization regime, the strength 
of the n = 2 signal is sensitive to the astrophysical details of the 
source, and especially to the black hole spin. In some cases, the 
n = 2 signal even vanishes entirely due to absorption. However, 
we find that at 345 GHz, the n = 2 signal returns to levels broadly 
consistent with previous estimates. This suggests that 345 GHz 
is an appropriate target for future space-VLBI observations of 
M87*, since the signal at that frequency is robust to the astro- 
physical parameters we vary. We stress that these results only 
mark the beginning of a systematic study of the reasonable pa- 
rameter space. In particular, our models exclude the possibilities 
of tilted disks, highly inclined observers, or jet-base emission. 


The paper is organized as follows. We present our accretion 
disk model in Sec. 2 and display images of some accretion flows 
in Sec. 3, before discussing the resulting visibility amplitudes in 
Sec. 4 and summarizing our conclusions in Sec. 5. 
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2. Accretion flow model 


In this paper, we consider accretion onto a black hole described 
by the Kerr metric in Boyer-Lindquist coordinates (t, r, 6, 6). The 
dimensionless spin parameter is labeled a. We use natural units 
such that G = c = 1, and radii are thus expressed in units of the 
black hole mass M. We always use p to denote the cylindrical 
radius p = rsin 6 (rather than a density) and z = rcos 6 to denote 
the cylindrical height above the equatorial plane 6 = 7/2. 


2.1. Disk geometry, physical quantities, and emission profile 


Our goal is to develop a very general and fully analytical model 
of geometrically thick disks, extending the one in Vincent et al. 
(2021). We restrict our attention to axisymmetric disks, so there 
is no dependence on ¢ in any of our physical quantities. 

We consider a population of thermal electrons that fills the 
spacetime outside the event horizon. The (fluid-frame) electron 
number density n.(p, z) is specified in the equatorial plane at the 


(cylindrical) horizon radius py = ry = 1+ V1 - a’, 


NesH = Ne(PH, Z = 0). (1) 
as is the electron temperature T.(p, z), 
Tex = Te(PH, Z = 0). (2) 


We define density and temperature profiles via a prescription 
similar to the one used in the analytical radiatively inefficient 
accretion flow (RIAF) models (see, e.g., Broderick et al. 2011): 


ro ae 2 
nino=neale) (spp) 
=] 
Tan) =Ten( =} : 


The dependence of the density profile on the cylindrical height z 
is thus a simple Gaussian with standard deviation s, = ap, where 
a is a free parameter that sets the opening angle of the accretion 
disk. Indeed, if we define the surface of the disk to be a cone of 
height s, above the equatorial plane, then a corresponds to the 
tangent of its opening angle. As a — 0, we recover equatorial 
models similar to those of GLM 20 and Paugnat et al. (2022), 
but here we keep a # 0 to explore the effects of disk thickness, 
which is expected to be significant for realistic RIAF flows (Yuan 
& Narayan 2014). 

The magnitude of the magnetic field, which is necessary to 
compute the synchrotron radiation, is prescribed by imposing a 
constant magnetization throughout the disk. That is, we fix the 
ratio of the magnetic field and particle energy densities 

2 
ae Br /4x (4) 


MpC>Ne 


(3a) 


(3b) 


where B is the magnetic field magnitude, m, the proton mass, 
and c the speed of light (which we have restored here for clarity). 
We include the effects of thermal synchrotron emission and 
absorption using the phenomenological expression derived by 
Leung et al. (2011). We present this expression and discuss its 
accuracy in App. A. As shown in App. B, the radial profile of 
thermal synchrotron emission in our model is approximately 


flr) & exp (<4) . (5) 


where ¢ is a model-dependent number that takes a typical value 
near 3 for our models. However, we emphasize that Eq. (5) is 
used for interpretation only; our computations are performed 
with the more precise expression (B.1) of Leung et al. (2011). 
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Fig. 1. Under the assumption of optically thin emitting matter concentrated very near the horizon, the range of reasonable appearances for models 
of accretion onto a Kerr black hole can be bracketed by two extreme toy models: equatorial, orbiting matter (brown), which produces a “wedding 
cake” structure (Gralla et al. 2019; Johnson et al. 2020), and spherical, infalling matter (purple), which produces a “shadow” (Falcke et al. 2000). 
Here, we show a 230 GHz intensity cut parallel to the spin axis of M87* (taken to be a rapidly spinning black hole with a = 0.94), including 
photons that orbit at most one full orbit around the black hole (up to n = 2 half-orbits), with each model normalized to its peak intensity. The 
brown and purple curves were ray traced from analytic models (see Sec. 3.4 for details), while the green curve is numerical data from Johnson 
et al. (2020) produced with the simulation pipeline described in Wong et al. (2022). The horizontal bars indicate the location of the Kerr critical 


curve and the apparent position of the equatorial event horizon. 


2.2. Dynamics 
2.2.1. Orbiting motion 


To compute images, we still need to specify the four-velocity of 
the emitting electrons in the accretion flow. We introduce a linear 
combination of an orbiting matter component and an infalling 
matter component that allows one to interpolate between the two 
extremes illustrated in Fig. 1. 

First, we define the orbiting component by an azimutal four- 
velocity field with vanishing radial and polar components. As 
discussed in App. C.1, it is nontrivial to find a prescription that 
is simple, everywhere smooth and becomes Keplerian at large 
distances from the black hole. We follow the prescription of Gold 
et al. (2020) and write the four-velocity 1-form as 


pe 


l+p’ 


uci’ dat = —urt* (— dt + £dd), f= (6) 


which is a special case of the general profile (C.16) described in 
App. C.1. This choice gives rise to a mild divergence Q ~ p~!/? 
at the poles, but is otherwise well-defined outside the horizon. 
The polar divergence is irrelevant for our models of thin and 
thick disks, as they have effectively no emission from the poles, 
and it has no noticeable effect even in the limit of fully spherical 
emission for infalling matter (Fig. 6 left). 


Unit-normalization of the four-velocity (6) requires 


' 1 
—ufite = (7) 


V—(g" — 2gE + 3907) 


Note that this circular motion is not geodesic. 


2.2.2. Infalling motion 


Next, we define the radially infalling matter component. Here, 
we simply assume that the motion is geodesic, with no azimuthal 
angular momentum and an asymptotically vanishing velocity. As 
reviewed in App. C, the resulting four-velocity takes the form 


u 
Ue rad 


oh O; + Ul yg Or + it O65 (8) 


Urad 


where 


weg = 8", yg =-VC1- 89", uh = 2%. 9) 


Note that the ¢ component is nonzero due to frame-dragging. 


2.2.3. Combined motion 


In this paper, we only consider purely circular motion or purely 
radial motion, but here we show that it is easy to combine them. 
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Symbol Value Property 
M 6.2 x 10° Mo compact object mass 
D 16.9 Mpc compact object distance 
a {0.01, 0.94} BH spin parameter 
@ = tan Op {0.1, 1} disk opening angle 
circ/thin {1.5,5} 
rad/thin {5,20} . 
Ne-H - - max density of electrons 
, circ/thick {0.7, 2} 
rad/thick {2,10} 
Ten 10! K max electron temperature 
o 0.01 magnetization 
i 160° inclination angle 
Vobs [230, 345] GHz observing frequency 
f 100 pas field of view 
NXN 20000 x 20000 image resolution 


Table 1. Parameters of our model. The maximum electron number density varies across models in order to ensure an observed flux of ~ 0.5 Jy 


at 230 GHz. It is expressed in units of 10° cm~3 


, with the first and second numbers referring to BH spins of a = 0.01 and a = 0.94, respectively. 


The models are labeled as follows: ‘circ’ means “circular rotation” (Sec. 2.2.1), ‘rad’ means “radial infall” (Sec. 2.2.2), ‘thin’ refers to a disk with 


opening-angle parameter a = 0.1, and ‘thick’ to a disk with a = 1. 


Indeed, one can linearly combine these motions to obtain the 
total four-velocity uv“ 0, = u' 0, + u' 8, + u® Oy. Following the 
notation of Pu et al. (2016), we introduce Q = u?/u' and write 


“= el — Br Utaa> 
Q= Ovire + qd — Bo Qraa = Qycire), 


where 0 < 6, < 1 and 0 < fy < | parametrize the superposi- 
tion of the circular and radial components (6) and (8). Note that 
only ut, is present in the first equation because our prescribed 
orbital four-velocity (6) has vanishing radial component. Unit- 
normalization fixes the time component of the four-velocity to 


wc fp bt sew? 
Bi + 281gQ + Bog Q?’ 
and its @ component is then given by u? = Qu’. 
Again, we will not examine this combined motion here, so 
the parameters 6, and £4 will not be discussed in the remainder 


of this paper—this section is simply meant to highlight that our 
model allows for very general flow dynamics. 


(10a) 
(10b) 


(11) 


2.3. Parameter choices 


In this paper, we are primarily interested in the effects of nonzero 
disk thickness, the motion of the emitting material, and the black 
hole spin. We thus consider extreme cases corresponding to a 
thin (@ = 0.1) or thick disk (@ = 1), purely azimuthal (Sec. 2.2.1) 
or purely radial (Sec. 2.2.2) motion, and low (a = 0.01) or 
high (a = 0.94) spin. The other source parameters are fixed 
to take likely—or at least reasonable—values for M87*. We 
choose an inclination of i = 160° and a mass/distance pair of 
M = 6.2 X 10° Mo and D = 16.9Mpc (EHT L5; EHT L6). For 
the magnetization o, we adopt a low value of 0.01, meaning that 
we consider a weakly magnetized disk.! The normalizations of 


' A commonly made distinction for black hole accretion flows is that 
between Standard and Normal Evolution (SANE) versus a Magnetically 
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the electron density and temperature are chosen such that the 
230 GHz flux is of the order of 0.5 Jy for all configurations, in 
accord with the analysis (and assumptions) of the 2017 EHT data 
(EHT L4; EHT L5). The model is illustrated in Fig. 2, and all the 
values of its parameters are given in Table 1. 


3. Image of the accretion flow 
3.1. Ray tracing and image orders 


We perform general-relativistic ray tracing in the Kerr spacetime 
with the model described in Sec. 2. We use the open-source code 
Gyoro (Vincent et al. 2011) to trace null geodesics backwards in 
time from a distant observer. The code integrates the radiative 
transfer equation along the null geodesics to evolve the specific 
intensity J, across the accretion flow. The output is an image, 
i.e., a map of J,. As discussed in App. D, the Gyoro integration 
parameters are chosen to ensure the highest possible numerical 
precision for the ray tracing and radiative transfer computations, 
while maintaining a reasonable computation time. 

We will frequently differentiate between image orders. The 
0'-order (n = 0) primary image is defined as the image produced 
by selecting the part of each null geodesic that extends from the 
observer screen to its first angular turning point (following the 
ray backwards into the geometry). The 1‘'-order (n = 1) image is 
produced by selecting the part of each null geodesic that extends 
between its first and second angular turning points, and likewise 
the 2"'-order (n = 2) image arises by retaining the contributions 
from the portion of each geodesic between its second and third 
angular turning points; we do not consider higher image orders. 
These image orders are illustrated in Fig. 2, and App. E describes 
the technical details of their computation. 


Arrested Disk (MAD). Prior investigations of high-order photon rings 
have focused on the MAD regime (Johnson et al. 2020; Chael et al. 
2021), and we are able to qualitatively reproduce these results using 
strongly magnetized (0 = 1) thin disks (@ = 0.1). Throughout this 
paper, we adopt a weaker magnetization 0 = 0.01 corresponding to the 
SANE regime (see, e.g., Fig. 1 of Porth et al. 2019). 
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Fig. 2. Density profiles of the ‘thick’ (a = 1, left) and ‘thin’ (@ = 0.1, right) disk models we consider. The plots display the poloidal (p, z) plane, 
with all azimuthal angles ¢ projected to one single point in the plane. The filled black area is the black hole event horizon. The solid black lines 
delineate the region in which bound photon orbits exist (the “photon shell”), with the circular-equatorial orbits (prograde and retrograde) indicated 
by black dots. The red color scale encodes the log-scale profile of the electron number density, with a floor set at 1% of the maximum density 
(this floor is only applied to this figure for readability; it is not applied in our model). The two solid red contours correspond to a density of 1% 
and 10% of the maximum. The dashed red lines enclose the points located at a height less than 3s, above the equatorial plane, where s, is the 
standard deviation of the Gaussian distribution controlling the electron number density above the equator [Eq. (3)]. In the right panel, we also 
show a high-order null geodesic in green, with blue letters marking the distant observer’s screen (S) and the two @ turning points (7, and T>) 
along the geodesic. Very sharp changes of direction appear at the 6 turning points. These are due to the 2D projection in the (p, z) plane of the 
three spatial dimensions of the null geodesic. The geodesic represented in 3D space would look perfectly smooth. The n = 0 part of the geodesic 
extends between S and 7), the n = 1 between 7; and T>, and the n = 2 between T> and the black hole. The yellow arrows highlight the portions 


of the null geodesic that are responsible for most of the n = 0, n = 1, and n = 2 emission. 


3.2. Resulting images 


The images of our models for spins a = 0.01 and a = 0.94 are 
presented in Figs. 3 and 4, respectively. All of them share the 
same qualitative appearance, and they all display the three main 
features of black hole images: 


e Acentral dark area, whose size depends on the astrophysical 
assumptions (e.g., Gralla et al. 2019; Chael et al. 2021). 

e A bright, narrow ring produced by strongly lensed photons 
that execute multiple orbits around the black hole on their 
way to the observer. This thin ring can be decomposed into a 
series of n = 1,2,... photon rings (often collectively referred 
to as “the photon ring”), each of which is a lensed image of 
increasingly higher order n of the accretion disk. These sub- 
rings, which are observable, exponentially converge to the 
theoretical critical curve, which is not (Bardeen 1973). The 
geometry of the critical curve is a pure function of the black 
hole mass and spin and of the observer inclination. However, 
the geometry and flux of observable photon rings still depend 
on the astrophysical assumptions (e.g., Paugnat et al. 2022). 

e A thick annular region of primary emission (produced by 
n = 0 photons travelling “straight” from the source to the ob- 
server, without orbiting around the black hole) that strongly 
depends on the astrophysical assumptions. In our images, it 
lies primarily within the thin photon ring. 


For each of the models in Figs. 3 and 4, we present the full 
image (left column) and images of the n = 1 andn = 2 rings only 


(second and third columns). The right column displays intensity 
along horizontal and vertical (i.e., perpendicular and parallel to 
the spin axis) cuts of the full image, zoomed in around the n = 2 
contribution, which is also displayed by itself. Clearly, the details 
of the n = 2 contribution are quite sensitive to the astrophysical 
assumptions. We now analyze the origin of these differences. 


3.3. Intensity contribution from the n = 2 ring 


Consider first the low-spin images, whose n = 2 contributions 
are shown in the rightmost column of Fig. 3. Two features are 
noteworthy: (i) the left and right peaks are of approximately 
equal intensity; (ii) the peaks in the models of radially infalling 
matter have lower intensity overall. Feature (i) is a consequence 
of the weak frame-dragging at low spin, which implies that the 
left and right geodesics follow a very similar path through the 
flow (Fig. 5 top panel). The suppressed intensity of the radial 
case relative to the circular case arises from a difference in their 
redshift factors that can be attributed to a complicated interplay 
of various competing effects, which we now explore in detail. 


To this end, we define a local increment of specific intensity 
loaded onto any given geodesic via radiative transfer. First, we 
introduce a local emissivity j, and absorptivity a,, as well as the 
Planck function B,, which depends only on the local tempera- 
ture and which equals (by virtue of Kirchhoff’s law) the ratio of 
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Fig. 3. Images and photon rings at 230GHz. Three leftmost columns: inner 50 jas of the low-spin (a = 0.01) brightness temperature maps. The 
two left columns share the same linear color scale, which goes up to a brightness temperature of 3.7 x 10!° K. The third column is in logarithmic 
scale, with the overall scale varying across panels to enable better visualization. The total specific flux of each image, as well as the maximum 
brightness temperature of each n = 2 image (in units of 10° K), are indicated in yellow font. The white dashed curve in the left column shows 
the primary image of the equatorial event horizon. Rightmost column: Horizontal (red) and vertical (blue) cuts of the full and n = 2 brightness 
temperature profiles, centered around the n = 2 peaks regions. The temperature ratios r between the two n = 2 peaks are provided. 


emission to absorption. We also need the optical depth 


(12) 


Ty.= { a, ds, 
optical path 


where the integral is taken over the portion of the ray connecting 
the source to the observer, so that exp(—t,) is the transmission. 
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The local intensity increment then reads 


ol, = iv [1 — exp (—a,, 5s)] x exp (-1,) x g? (13a) 

= B,[1- exp (-a, 6s)] x exp(-ty) x g? (13b) 
ea OTe eee Ne 
self-absorbed local emission transmission __ redshift 

~ By exp(-Ty) 8°, (13c) 


where 6s is the increment of length along the ray as measured by 
the emitter, and g denotes the redshift factor 
_ Uobs * Pobs 


Uem * Pem 


(14) 
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Fig. 4. Same as Fig. 3 for high spin (a = 0.94). The radial-thin and radial-thick horizontal cuts do not show the left peak of the n = 2 ring because 
it is too small to be visible (the associated temperature ratios are r = 10* and r = 108). The right column has a different scale compared to Fig. 3. 


Here, Uobs/em Aenotes the four-velocity of the observer/emitter 
and Pobs/em the photon momentum at the observer/emitter. All 
the quantities appearing in Eqs. (13) depend on the local value 
of the emission frequency, which itself depends on the observed 
redshift, which in turn depends on the flow dynamics. 

To illustrate this, we show in the top row of Fig. 5 some of 
the geodesics that most contribute to the n = 2 peaks visible in 
the rightmost column of Fig. 3. For the circular/thin model, these 
light rays collect much more intensity than the analogous light 
rays for the radial/thin model, whose contributions are greatly 
suppressed by the redshift factor. This explains point (11) above. 

To summarize, the heights of the n = 2 peaks appearing in 
the rightmost columns of Figs. 3 and 4 are controlled by three 
factors that contribute to the RHS of Eqs. (13), namely: 


1. the emission, which in our model depends mostly on radius, 
as itis highly concentrated near the equatorial plane (near the 


n = 2 equatorial crossing of each light ray); however, it also 
depends on frequency and therefore redshift (see 3. below); 

2. the absorption and transmission, which behave similarly to 
the emission, with the important difference that they are non- 
local, as they depend on the full history of the geodesic as it 
travels from the n = 2 emission site to the observer; 

3. the redshift, which has a complicated dependence on radius 
(the gravitational redshift is linked to the metric and blows 
up at the horizon), as well as on the relative orientation of 
the emitted photon and emitting particle velocity. 


The interplay of these three ingredients strongly depends on the 
observing frequency, on the dynamics of the flow, on the lo- 
cal physical conditions at the emission site (density, magnetic 
field, temperature), as well as on the physical conditions along 
the complete path followed by the photons. It is therefore not a 
simple task to predict the properties of the n = 2 contribution, 
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Vem=690 GHz 


~ Vem=275 GHz, em=2.5e-7 
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Vem=275 GHz, em=2.6e-7 
tr=0.0025, g?=0.57, 61,=3.7e-10 


Vem=2.9 THz 
em=3.8e-5 
tr=1e-6 
g?=5e-3 
61,=1.9e-13 


Fig. 5. Origin of n = 2 emission. We show geodesics corresponding to the maximum intensity of the n = 2 ring of each model along a horizontal 
cut, on the left (blue) and right (red) sides of the image (see inset of the upper-left panel). The accretion flow and spin parameters are specified in 
the upper-left corner of each panel. The blue and red ellipses encircle the region emitting most of the n = 2 photons loaded onto each geodesic. 
The blue geodesic is abruptly cut in the lower-right panel because the medium becomes optically too thick (defined in the code as a transmission 
smaller than 10~°). For each ellipse, the local values of the emitted frequency (Vem), self-absorbed emission (em, in cgs units), transmission (tr), 
and redshift factor (g?) are provided, as well as the resulting increment of specific intensity 5/, (see Eq. (13) for the definition of these quantities). 


and simulations are certainly the only way to investigate them. 
Figure 5 provides the typical values of these three quantities in 
the region of n = 2 emission for each of our models. 


At high spin, the n = 2 peaks (shown in the rightmost column 
of Fig. 4) are very different from their nonrotating counterparts. 
They also display two striking features: (i) in stark contrast with 
the nonrotating case, the left and right peaks in the horizontal 
cuts (perpendicular to the spin vector, shown in red) have very 
different heights, indicating a large brightness asymmetry in the 
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image; (ii) this brightness asymmetry ratio (the intensity ratio of 
the two horizontal peaks) is extremely high for the radial models. 

These features can be understood by examining the bottom 
row of Fig. 5, which shows that the left and right geodesics in 
the high-spin case probe very different parts of the inner flow, 
resulting in a different collected intensity. This is closely related 
to the fact that as black hole spin increases, the critical curve 
is both distorted (from a circle) and displaced (relative to the 
primary radiation, see, e.g., Chan et al. 2013, Fig. 5 left panel): 
since the n = 2 photon ring shares these properties, the left and 
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Fig. 6. Intensity cuts along the directions perpendicular to spin (red) and parallel to spin (blue) for spin a = 0.01 for the various thick-disk models 


used in this article (four central panels), as well as for a spherical model 


obtained from the thick-disk model in the limit of large disk thickness 


(top left panel). The size of the critical curve (thin line) and of the primary (n = 0) image of the equatorial horizon (thick line) are shown in green. 


right geodesics visit different parts of the flow, which in turn 
leads to the large horizontal (perpendicular-to-spin) brightness 
asymmetry of the n = 2 ring. This explains point (i) above. 


Note that for the vertical cuts (parallel to the spin vector, 
shown in blue in the rightmost column of Fig. 4), the brightness 
asymmetry is much milder. This makes sense because at high 
spin, the distortion and displacement of the critical curve and 
n = 2 ring primarily act in the direction orthogonal to the spin. 


The extreme horizontal brightness asymmetry obtained in 
models of radial infall can be attributed to a strong absorption of 
the n = 2 radiation along one of the geodesics, as well as to an 
extreme redshift effect due to the fact that the n = 2 emission is 
loaded onto the left geodesic (in blue in Fig. 5 lower-right panel) 
very close to the event horizon. This explains point (ii) above. 


3.4. Intensity cuts 


We conclude this section by revisiting the framework proposed 
in the introduction, in which realistic models are regarded as 
lying on a continuum between the shadow and wedding cake 
extremes. To this end, here we discuss intensity cuts of all the 
models analyzed in this paper, starting with the extreme models 
(spherical and equatorial) and then moving on to our thick-disk 
models. 


3.4.1. Extreme models: spherical & equatorial 


Figure | of the introduction presents the intensity cuts of two 
analytical models that are not based on the thick-disk model that 
we develop here. These two models are defined as follows: 


e The spherical infall model consists of a spherically symmet- 
ric distribution of emitting matter located everywhere outside 
the event horizon. This matter is falling into the black hole 
with four-velocity defined in Sec. 2.2.2. The number den- 
sity and temperature of the emitting fluid are chosen at the 
event horizon and evolve with radius following the power 
laws defined in Eqs. (3) (i.e., ne(r) r?, T(r) « r7!). The 
magnetic field is prescribed by demanding that the magneti- 
zation o defined in Eq. (4) remain constant, at the same value 
as in our thick-disk model, @ = 0.01. We consider thermal 
synchrotron emission and ignore self-absorption. 

The equatorial orbiting model is that of GLM 20 with an ad- 
hoc emission profile (their emission model 1). There is also 
no absorption considered in this model. 


3.4.2. Canonical thick-disk models 


We now discuss the intensity cuts of our four canonical thick- 
disk models, which are displayed in Fig. 6. We focus on the case 
of low spin, as the high-spin case is qualitatively similar. 

The “spherical limit” in the top-left panel is the a — oo 
limit of our disk model with radially infalling matter. This large- 
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thickness limit of our thick-disk model exactly coincides with 
the spherical infall model discussed in the previous subsection. 
The only difference is that in the upper-left panel of Fig. 6, 
we take synchrotron self-absorption into account. This spheri- 
cal limit is similar to the models used by Falcke et al. (2000) and 
Narayan et al. (2019), except that we self-consistently include 
absorption in our model of synchrotron emission. 

The a — 0 limit of our thick-disk is numerically tricky, and 
the emission profile of the analytical equatorial model depicted 
in Fig. | differs from that obtained with our thick-disk models. 
For this reason, we do not show an “equatorial limit” intensity 
cut in Fig. 6. 

Figure 6 shows a continuum of models spanning a large pa- 
rameter space, from near-equatorial to spherical emission, and 
from circularly orbiting matter to radially infalling matter. We 
highlight a few important observations: 


e The radial spherical model (top left) shows a prominent pho- 
ton ring surrounding a large shadow that essentially occu- 
pies the full interior of the critical curve. The existence of 
the shadow is due to the strong redshift effect on the radia- 
tion emitted by matter that falls towards the black hole.” The 
spherical model is the only case where the intensity profile 
does not decompose into discrete layers (n = 0, 1, 2). 

Our thick-disk models produce a wide range of profiles that 
interpolate between the spherical limit and the equatorial 
model depicted in Fig. 1. Like the equatorial model, these 
models all produce intensity profiles that decompose into 
distinct n = 0 and n = | contributions. The n = 2 photon ring 
is also clearly present for all the circular models; in the radial 
models, it is also present but (as discussed above) weaker and 
less easily discerned. The contrast between the photon rings 
and the n = 0 emission is stronger for the radial models be- 
cause their n = 0 contribution is strongly suppressed (by the 
same redshift effect that produces a shadow in the spherical 
model). As a consequence, the central brightness depression 
is larger in the radial models than in the circular models; for 
the latter, as for the equatorial model, the central dark region 
is restricted to the interior of the apparent equatorial horizon. 


4. Visibility signatures 


One of the motivations for this paper is the prospect of an exper- 
imental detection of the n = 2 photon ring via its long-baseline 
“ringing” in the Fourier plane, which might be observable with 
future space-based VLBI missions (Johnson et al. 2020; Pesce 
et al. 2019; Gralla 2020; Gralla & Lupsasca 2020c; GLM 20). To 
check whether such a feature would be observable in our models, 
we compute the visibility along cuts in the Fourier plane in the 
directions parallel and perpendicular to the spin axis. Following 
GLM 20, we leverage the projection-slice theorem to simplify 
the calculation. For a given orientation, we determine the line in- 
tegrals on all lines perpendicular to the orientation (a 1D cut of 
the Radon transform of the image), and we Fourier transform this 
1D function to find the desired cut of the 2D Fourier transform, 
which is the radio visibility. The magnitude of this 1D complex 
visibility is the visibility amplitude, which we will simply refer 
to as the visibility. In the same way that we extracted from the 
full image its n = 1 and n = 2 components, we will compute the 


> Null geodesics arriving inside the critical curve have no radial turning 
points so the emitted radiation is directed towards increasing radii, while 
the emitter is heading towards decreasing radii, giving rise to strong 
redshift suppression of the observed intensity via the “headlight” effect. 
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full Fourier transform and its n = | and n = 2 components. Our 
Fourier transform conventions agree with GLM 20. 

Table 2 shows the main features of the visibilities associated 
to the various models considered in this study. In particular, it 
provides the visibility amplitude level at 40 and 100. GA, as well 
as the baseline thresholds (and corresponding visibility levels) 
past which the n = | and n = 2 rings start to dominate the signal. 
These thresholds were called b; and b2 in Paugnat et al. (2022). 

The following criterion is used to define when the n" ring 
starts to dominate the visibility. The full image decomposes into 
layers indexed by n, /,(x) = Ip(x) +/,(k) + h(x) +..., and the full 
visibility V(u) = Vo(u) + V; (a) + V2(u) +... does too, by linearity 
of the Fourier transform. We consider a sliding baseline window 
u ~ Uy Of fixed width 10GA, which is approximately twice the 
photon ring diameter. At each angle yg in the baseline plane, we 
compute over this window u ~ u,, the mean visibility (|V(u, ¢))),, 
of the full image /,(x) and the mean visibility (|V,,(u, y)|),, of the 
n' ring image J,,(x) only. We define the baseline threshold b,(¢) 
past which the n" ring dominates the visibility to be the shortest 
baseline window u,, = b,(y) such that the percentage difference 
between these two means dips below a given threshold p: 


Vi Ow — (VU, Py n ring dominates 


(VU Pw 


(15) 
onu z bily) = Uy. 


We will generally take p = 5%, unless otherwise stated. 

We again focus our comments on the n = 2 ring signature. 
In all low-spin models, the n = 2 ring clearly dominates the 
signal on sufficiently long baselines u 2 b2, with the threshold 
by ranging from ~ 600 — 1400.GA. The corresponding visibility 
amplitude displays clear oscillations at levels varying between 
~ 10 — 50yJy (see Table 2 and Fig. 7). Thus, these low-spin 
models match the idealized models considered in GLM 20 and 
Paugnat et al. (2022) reasonably well. 

At high spin, by contrast, it is only on the baseline oriented 
parallel to the spin axis (y = 90°) that the n = 2 ring clearly dom- 
inates the signal (see the four lower panels of Fig. 8). However, 
even at this orientation, the models of radially infalling matter 
produce a very weak visibility amplitude in the n = 2-dominated 
regime ~ 0.01 —0.1 wJy. This is directly related to the low height 
of the n = 2 peaks in the rightmost column of Fig. 4: the blue 
profiles for the two models of radial infall display a contrast with 
the n = | radiation of ~ 100 — 1000. 

On baselines perpendicular to the spin axis (y = 0°), all the 
high-spin models produce a very weak n = 2 signal, with essen- 
tially no oscillation at long baselines (see four upper panels of 
Fig. 8). This is the reason why some numbers are missing in Ta- 
ble 2: these entries correspond to cases in which the n = 2 oscil- 
lation is vanishingly small. To understand this effect, recall that 
the oscillations in visibility amplitude are caused by “interfer- 
ence” between opposite peaks in the image, with the amplitude 
of the oscillation set by the smaller of the two (Eq. (7) in Gralla 
2020). Given that all high-spin, spin-perpendicular cuts display 
a large brightness asymmetry between the two n = 2 peaks (see 
the red profiles of the rightmost column in Fig. 4), the resulting 
visibility cuts exhibit essentially no oscillations associated with 
the n = 2 ring. This asymmetry between the n = 2 peaks is itself 
due to the interplay between absorption and redshift effects on 
the various parts of the n = 2 ring (see Sec. 3.3), as illustrated by 
the geodesics plotted in the bottom row of Fig. 5. 

To summarize what we have discussed so far, the n = 2 peaks 
in the intensity cuts displayed in the rightmost column of Figs 3 
and 4 have two features of fundamental importance: (i) the ra- 
tio of flux in the n = 2 image /)(x) relative to the total flux in 
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Model 40GA 100GA n=1 n=2 
(150Ga, (1400 Ga, 
_ : ode 9.2 mJy 2.4mJy I mJy) 25 ply) 
@=00)cucalartin: | 49sniy 2.7 mJy (140Ga, (1300Ga, 
2 mJy) 20 py) 
(54Ga, (800 Ga, 
4 = 0.01 radial thin on 14.4 mJy saan 2.2 mJy sea 10 mJy) ewer 10 pJy) 
7.5 mJy) 5.6 wy) 
(54 Ga, (650GaA, 
_ e , 8.6 mJy 2.4mJy 6.8 mJy) 55 py) 
Pe OO Lorena || eae, 2.8 mJy (44.Ga, (650GaA, 
7.71 mJy) 43.2 uy) 
(35 GA, (660 Ga, 
4 = 001 tadial thick eeu 12.2 mJy pie 2.2 mJy en 14.4 mJy) , 23 py) 
12.8 mJy) 14.9 wy) 
(i2Ga, (700 Ga, 
ee ce een oe 12.0 mJy one 2.0 mJy aie. 1.8 mJy) aan. 57 pty) 
N/A) N/A) 
(44.2 GA, (896 GA, 
psirOA adhe an 10.4 mJy ane 1.3 mJy eae 9.7 mJy) sine 0.25 py) 
5.5 mJy) N/A) 
(54GA, (328 Ga, 
_ : : 10.7 mJy 1.8 mJy 6 mJy) 180 py) 
BE CUI MEI| e aey 0.7 mJy (63.8GA, (N/A, 
3.5 mJy) N/A) 
(34.4 GA, (847 GA, 
_ 343 . 7AmiJy 0.6 mJy 10.0 mJy) 0.01 py) 
a= 0.94 radial thick | 6 9 mJy 0.5 mJy G4.4Ga, (N/A, 
7.3 mJy) N/A) 


Table 2. Key features of the visibility amplitudes of the various models considered at v,,, = 230GHz: for each model, we list the visibility 
amplitude |V(u, y)| at u = 40 GA and 100 Ga, as well as the baseline thresholds b;(y) and b2(y) (and the corresponding visibility amplitudes) past 
which the n = 1 and n = 2 contributions start to dominate the signal [Eq. (15)]. Each box contains two sets of numbers separated by a diagonal: 
these correspond to a baseline polar angle y = 0° (perpendicular to spin, lower left) or y = 90° (aligned with spin, upper right). The percent level 
p of Eq. (15) is 5% everywhere, except for the numbers in red, where it is set to 10% (in these cases, the ring never dominates according to the 5% 
criterion). The ‘N/A’ in the n = 2 column for the high-spin cases at y = 0° mean that the n = 2 signal never dominates even for p = 10%. 


the full image /,(x) controls the average strength of the visibility 
in the n = 2-dominated regime; (ii) the brightness asymmetry 
(intensity ratio) between the left and right (or top and bottom) 
peaks in the n = 2 image cuts controls the amplitude of the vis- 
ibility oscillations in the n = 2-dominated regime. These em- 
pirical observations are in perfect agreement with the analytical 
study of Gralla (2020). They allow one to guess whether a model 
will produce a strong n = 2 signal by simply computing a high- 
resolution intensity cut along a given orientation in the image 
plane, which is very cheap in terms of computational resources. 

The fact that the n = 2 signal is absent from the spin- 
perpendicular visibility in the high-spin models is problematic 
from an experimental perspective. However, since absorption 
is the main culprit in suppressing the n = 2 contribution, one 
might expect the signal to reappear at higher observation fre- 
quencies, where absorption is lower. While we have so far re- 
stricted our attention to the present-EHT observation frequency 
Vobs = 230 GHz, we can also consider a higher frequency vops = 
345 GHz of planned future VLBI observations (Doeleman et al. 
2019; EHT L2). Table 3 provides the same information as Ta- 
ble 2, for the high-spin models observed at vop; = 345 GHz. This 
higher frequency contains a clear n = 2 signal in all cases, and 


the associated visibility level is in the regime ~ 100 Jy that 
is considered potentially observable by GLM 20. The stronger 
n = 2 signal at higher frequency is also readily seen in the im- 
age cross-sections. For example, Fig. 9 displays intensity cuts for 
the circular/thin model at these two frequencies, and the higher- 
order peaks are clearly much sharper at 345 GHz. 


5. Conclusion 


In this article, we studied the images and visibility amplitude 
profiles of a variety of models of geometrically thin and thick 
disks, including the effects of absorption. We found that the size 
of the central brightness depression depends strongly on astro- 
physical assumptions, and that the photon ring is generically dis- 
crete rather than smooth, supporting the “wedding cake” heuris- 
tic over the “black hole shadow” paradigm. We also studied the 
long-baseline visibility amplitude, focusing on the signature of 
the n = 2 photon ring of M87*. Although the n = 2 signal is al- 
ways clear at low spin at 230 GHz, it can disappear at high spin 
due to a combination of frame-dragging and absorption along 
the line of sight. However, at 345 GHz, we find a clear n = 2 sig- 
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Fig. 7. Visibility amplitude profiles for all low-spin (a = 0.01) models. The full profile |V(u, y)| is shown in black, the n = 1 profile |V\(u, y)| in 
green, and the n = 2 profile |V2(u, y)| in red. The four upper panels correspond to an orientation y = 0° (perpendicular to the spin axis, -L), and the 


four lower panels to an orientation y = 90° (parallel to the spin axis, ||). The accretion flow model is specified in the top-right corner of each panel. 
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Fig. 8. Same as Fig. 7 for the high-spin (a = 0.94) models. In all models, the n = 2 signal is very weak at y = 0° (orientation perpendicular to the 
spin axis) due to the strong frame-dragging effect at high spin (see text for details). Note the different vertical scale for the two lower right panels. 
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Table 3. Same as Table 2, but at y9,, = 345 GHz. Only high-spin models are considered. 


nal even at high spin, because the absorption is less prominent at 
this higher frequency. 


We thus find that 345 GHz is a promising target for future 
space- VLBI observations of the photon ring of M87*. However, 
we emphasize that this conclusion applies only for the particu- 
lar astrophysical conditions that we have considered herein. Al- 
ternative reasonable choices of inclination, total compact flux, 
magnetization, density profile, etc., may very well lead to a dif- 
ferent conclusion. We also cannot discount the possibility that 
the source geometry is qualitatively different, consisting for in- 
stance of a tilted disk (e.g., White et al. 2020), a jet-dominated 
profile (e.g., Kawashima et al. 2021), or emission from current 
sheets (e.g., Crinquand et al. 2022). Finally, we should keep in 
mind that the source itself may be variable in ways that we do 
not understand and cannot predict. For example, the total com- 
pact source flux density appears to have varied by a factor of 2 
across a decade of VLBI observations, which would influence 
the density scale and hence the system optical depth (Wielgus 
et al. 2020). It is therefore of paramount importance to extend 
our analysis to a larger astrophysical parameter space covering 
all reasonable possibilities for the source. 
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Fig. 9. Intensity cuts along directions perpendicular (red) and parallel 
(blue) to spin for the circular/thin model at high spin, for two observing 
frequencies: 230 GHz and 345 GHz. The emission region is more opti- 
cally thin at 345 GHz, so the higher-order peaks are more pronounced. 
In particular, the n = 2 left peak of the spin-perpendicular case is very 
sharp at 345 GHz, and fully absorbed at 230 GHz (see green ellipses). 
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Appendix A: Relativistic thermal synchrotron 
emission and absorption 
The synchrotron power per unit frequency emitted by a single 


ultrarelativistic electron is well-known (see Rybicki & Lightman 
1979, Sec. 6.2) to be 


dE eBsind _( v 
= V3 F : A.l 
dt dv mc? Veit acl) 
where 
3 4 ; eB 
Verit = a” Veyclo SIN 0, Veyclo = Sane (A.2) 


Here, y denotes the Lorentz factor of the electron, 6 the angle 
between the direction of emission and the magnetic field, and 


0 2/37 (2) +1/3 
Foyax [ Keptodu= | r(3)x'? forx <1, (A.3) 


Vue * for x > 1, 


with Ks/3(u) a modified Bessel function of the second kind. Note 
that in these formulas, we have identified the pitch angle between 
the magnetic field and the direction of the electron’s velocity 
with the emission angle 6, which is a valid approximation for 
ultrarelativistic electrons (Rybicki & Lightman 1979). Also note 
that the derivation of Eq. (A.1) is essentially built on the strong 
beaming effect of the gyrating ultrarelativistic electrons, which 
is thus taken into account from the very start in this treatment. 
The associated thermal synchrotron emissivity is 


dE 3B si . dng 

iz _ N3e =| F v Feiss (AA) 

dtdvdVdQ 4x mc? 0 Verit} dy 
where 

x 1/2 

dne WPT) ( Y 

= exp|— (A.5) 
dy e. K> (0.') : ©. 


is the relativistic thermal (Maxwell-Jiittner) distribution, K2(u) is 
a modified Bessel function of the second kind, ©, = kgT./ mc? is 
the dimensionless electron temperature, and the factor of (4x)! 
implicitly assumes an isotropic distribution of the electrons’ mo- 
menta. Note that the integral over y in Eq. (A.4) is evaluated over 
the range [0, co), even though physically, y > 1. Extending the 
range of integration to include the interval [0, 1] simplifies the 
computation of the integral, while maintaining a good accuracy 
(as we will see below). Using the asymptotic expansions of F(x) 
for x < 1 and x > 1, the emissivity (A.4) may be analytically 
expressed as (see, e.g., Leung et al. 2011) 


4/3 
24/37 nee2V,y 


5 e ax for X « 1, 
: c 
y~ A. 
J V2re? Vs —x1/3 ( 6) 
Ne “eaae- for X > 1, 
Cc e 
where 
Vv 2 Di 26 
X=-—, Vs = 9 cycloOs Sin 6. (A.7) 
Vs 


Leung et al. (2011) also provide a fitting function that bridges 
these two asymptotic regimes: 


V2ne"y, 
*3cK> (@;') 


approx 


i =n (x? ae guney 5), exe (A.8) 
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To summarize, Eq. (A.8) is derived under the following list 
of approximations: we assume the electrons are ultrarelativistic 
and their population isotropic, we slightly extend the integration 
bounds on y, and we match the asymptotic expansions of F(x). 

In our ray tracing code, we use the formula (A.8) for the 
emissivity, averaged over the emission direction 6: 


a) = i ie dQ = 1 a i hee sin dd, 
7 2 Jo 


where dQ = sin @d6d¢ and the emissivity is independent of the 
azimuthal angle @. 

The exact expression for the emissivity, free from the above 
approximations, is (see Leung et al. 2011 for details) 


(A.9) 


2ne?v YM) 7 dn 
ot 7 € A.10 
or 2 Ts | plas (A.10) 
where n_ = 
Veyel Veyelo \2 a) 
n—— + |cos 6 (n=) — sin’ @ 
y(n) = ; (A.11) 


sin? 6 
where n is the integer that is being summed over in Eq. (A. 10). 


The dimensionless function J,,(y) is given in terms of the Bessel 
function of the first kind J,,(z) and its derivative J/(z) by 


Tn) = MIP + [NI@), (A. 12a) 
(4) = 
go = YF sin Osinag, (A.12b) 
sin 0 Veyelo 
1 cyclo 
N =Bsina, COS Qo = =o : (A.12c) 
Bcos@ y v 


The sum over 7 in Eq. (A.10) accounts for the helical motion of 
electrons along magnetic field lines, which makes each electron 
emit at every multiple of its Doppler-shifted gyrofrequency. 

The accuracy of the approximation (A.8) can be assessed by 
comparing it with the exact expression (A.10). This is done in 
Figs. 10 and 11 of Leung et al. (2011), where the relative error 
between these two quantities at 6 = 30° and 80° is displayed as 
a function of @, and v/Vcyclo. This error is typically of the or- 
der of a few percent for ©, € [1,10] and v/Veyc1o € [107, 105}, 
which includes the typical ranges for these parameters in the in- 
ner region of the accretion flow surrounding M87*. Indeed, the 
temperature of the source is T, ~ 10'' K, so @, ~ 15 (and the ac- 
curacy gets better with increasing ©), while the magnetic field is 
B = 10G, which implies a cyclotron frequency Veycio © 10’ Hz, 
so that v/Veyclo © 10*. On the other hand, the error typically ex- 
ceeds 10% when ©. < 1, which means that the approximation 
(A.8) does not accurately model the emission from the outer re- 
gion of the accretion flow. Nevertheless, this is not an issue for 
us because these contributions are negligible compared to the 
emission from the inner region of the flow. 

The thermal absorption coefficient is simply related to the 
emissivity by Kirchhoff’s law a, = j,/B,, where B, is the Planck 
function, so the approximations for the absorption are exactly the 
same as for the emission. 


Appendix B: Synchrotron emission profile for M87* 
Appendix B.1: Thermal synchrotron 


Since ©. ~ 15 for M87*, we may assume that ©, > | for this 
source. In that case, the approximation (A.8) for the emissivity 
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Fig. B.1. Ratio of (X!/? + 2!!/!?X'/6)? to X [or equivalently, of the emis- 
sivity (B.1) to its approximation (B.3)] as a function of magnetic field 
strength B and dimensionless temperature ©,, fixing sind = 1 and 
v = 230 GHz (and neglecting redshift effects). The solid white curves 
are level sets of ratios of 1.5, 2 and 3. The dashed red curves are level 
sets of the X parameter corresponding to X = 10 and X = 100. 


; : : 0 : bg died 
can be further simplified using K>(x) “x 2x2. This substitution 
results in the further approximation 


V 2767 Veycto sin 6 


2 _yi/3 
1/2 11/12 y1/6 x 
ae (x apuey ye , 


Ty SMe (B.1) 
which is accurate to better than 10% provided that ©, > 1.5, and 
to better than 1% provided ©, > 5. 


Plugging typical values for M87* yields an estimate for X of 


y ~ 1600 ( vfHzI\/_10_\/ 10)’ 
~ sind \ 102 )\B[G]]\@.} ’ 
where the frequency v and magnetic field strength B should be 
expressed in Hertz and Gauss, respectively. The factor of 1600 
incorporates all the constant terms in the definition (A.7) of X. 
Since this expression satisfies X >> 1 for all typical values 
of M87* parameters, we can further simplify the approximation 
(B.1) for the emissivity by taking its X — oo limit. In the regime 
X > 1, we may replace (X!/* + 2!!/!2X1/6)? by X to find 


(B.2) 


. O6,.X>1 


V2me7 Veyclo INO _yiss 
Jv Ne ———.—_ Xe 


B. 
27c , — 


which is a slight underestimation because (X!/? + 2!!/1?X1/6)? is 
always bigger than X. Their precise ratio is plotted as a function 
of B and ©, in Fig. B.1 and is always of order ~ 1 — 10 for typical 
parameter values for M87*. 

Next, since we are only interested in the overall profile of j,, 
we may drop all constant prefactors and write 


. 1/3 
(se) | 
Bez 


where C = (92mc/e)'/? is a constant. 


(B.4) 


Since we are assuming that the density and temperature fall 
off as r~-? and r~', respectively, we may immediately conclude 
that the radial dependence of the emissivity is simply given by 


jr) & exp [-X(r)"], (B.5) 
provided that we neglect the radial dependence of the frequency, 
which is affected by a position-dependent redshift; this is only a 
reasonable assumption provided the emission is not too close to 
the horizon r = r,, where the redshift diverges. 

Returning to the estimate (B.2) for X, and using the assumed 
magnetic field fall-off B « r} [Eq. (4)], we may now write 


XxX & 


5 3 
3.7 x 10 ( r (B.6) 


2 
B innerO 


e;inner sin 6 Pinner 
where we have set vy = 230 GHz and once again neglected the 
frequency redshift. Here, inner 1s the innermost radius of the disk, 
while Binner and ©¢ inner respectively denote the innermost values 


of the magnetic field and dimensionless temperature. Hence, 


1/3 
_ (B7) 


To summarize, Eq. (B.7) is derived under the approximations 
listed below Eq. (A.8), together with the following additional as- 
sumptions: we assume that ©, >> | and X > | (which puts us 
between the red contours of Fig. B.2), and we prescribe power- 
law fall-offs for the density, temperature, and magnetic field as 
described in our model (n, « r-?, @ « r7!, B o r7!). Should 
the various power-law indices differ from this choice, one could 
derive a corresponding correction to Eq. (B.7) starting from the 
more general expression (B.4). We reiterate that this approxima- 
tion underestimates the emissivity by a factor of a few, and that 
its quality increases with X (Fig. B.1). 

The thermal-synchrotron emission profile for typical M87* 
parameters thus decays exponentially in the radius, with a slope 
¢ depending on the exact conditions. Table B.1 lists the values of 
¢ for our geometrically thin models, and Fig. B.2 displays ¢ as a 
function of the magnetic field strength and temperature. 

For thermal synchrotron emission, the absorptivity may be 
obtained from Kirchhoff’s law as 


: r 3.7 x 10° 
Iv) oc eXp -- | > g _ Fer 


e;inner sin 6 


Aa = jy 
"BY. 


(B.8) 


where B,(T.) is the Planck function at the local electron temper- 
ature T,. At 230 GHz, photons are deep in the Rayleigh-Jeans 
regime where hy « kpgT-, so that this relation simplifies to 


2 


c 
y R —— J. B.9 
Oy kT ev?” po 
Assuming as usual a temperature profile T « r~', we obtain 
a@,(r) « exp -< | : (B.10) 
Tinner inner 


The highly lensed emission originates from small radii r 2 inner, 
for which the exponential profiles (B.7) and (B.10) provide good 
approximations to the emissivity and absorptivity, respectively. 
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Fig. B.2. Values of the parameter ¢ [Eq. (B.7)] as a function of the 
magnetic field strength B and dimensionless temperature ©,, fixing 
sin @ = 1. Our models always stay in the part of the plane located above 
the dashed horizontal red line (©. > 5, where the approximation (B.1) 
is good) and below the dashed red contour (X > 10, where the approx- 
imation (B.3) is excellent). The solid white contours are level sets of 


é. 


Model ra 
a=0.01 circ/thin | 4.0 
a=0.01 rad/thin | 3.4 
a= 0.94 circ/thin | 3.4 
a=0.94 rad/thin | 2.7 


Table B.1. Values of the parameter Z in a few of our models, for which 
the emissivity profile is approximately exp(—¢r/Trinner), provided that r 
is far enough from the event horizon and redshift effects are neglected. 


Appendix B.2: Power-law synchrotron 


Although we only study thermal-synchrotron emission in this 
paper, it is interesting to also consider a nonthermal population 
of electrons (accelerated by magnetic reconnection, for instance) 
with synchrotron emission following a power law with index p. 
Still following Leung et al. (2011), the emissivity is then 


-PL € Vey 3P/2(_p = 1) sin @ 
Jy = Ne I-p _ ,|-p 
c 2(p + 1) (ae ~ Venu | 


er ( 321) (32419) (_» Soe 
12 1) eae 


where Ymax/min are the maximum/minimum values of the Lorentz 
factor for electrons in the distribution and I(x) the gamma func- 
tion. Following the same steps as in the previous section leads 
to 


cara : 


(B.11) 


; (B.12) 


Tinner 


i 2 


i.e., a power-law emissivity profile depending only on the power- 
law index of the electron distribution, and not on the physical 
conditions of the flow. For the particular value p = 3, jP/ « r-?. 
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The power-law absorptivity a?" is also provided by Leung 
et al. (2011) and also follows a power law that depends only on 
the power-law index of the electron distribution. For p = 3, it 
behaves as a?! « r°/*, Thus, it appears that the radial profile 
of the emissivity and absorptivity can vary significantly with the 
choice of electron distribution, and this could in turn have signif- 
icant impact on the n = 2 ring signature (see Sec. 3.2). It would 
therefore be useful to also investigate this signature in models 
with nonthermal emission. 


Appendix C: Geometrically thick accretion flows 
Appendix C.1: Orbiting motion 


In this section, we restrict ourselves to the Schwarzschild metric 
for simplicity. We want to prescribe a four-velocity for circularly 
orbiting matter in a thick disk. A circular four-velocity u“ and its 
associated 1-form uw, take the general form 


Ws Oy = ul (A, + Qs), (C. 1a) 
us" dx" = —u, (— dt + €d¢), (C.1b) 
where Q = u?/u' denotes the angular velocity and € = —ug/u; 


the specific angular momentum of the flow. 


Appendix C.1.1: Equatorial Keplerian motion 


Timelike circular geodesic motion in the equatorial plane (i.e., 
Keplerian motion) is well-known (Bardeen 1973). The four- 
velocity takes the form (C.1) with 


r-2 r 
Wir) 


It is obvious from these relations that Keplerian motion is not 
allowed inside the photon shell at r = 3. Moreover, Keplerian 
motion is only stable outside the innermost stable circular orbit 
(ISCO) located at r7sco = 6. 

Below the ISCO, Cunningham (1975) prescribes that the 
geodesic constants of motion —u,; and uy (the energy and spin 
angular momentum) keep their ISCO value, which ensures that 
the flow remains continuous across the ISCO as it spirals into the 
horizon. This prescription results in the four-velocity 


0 = r3??, 


(C.2) 


—Uuy; = 


uy, dx = uP? dt + u,dr + up dé, (C.3) 


where the (t, #) components ee = 


t ug(r = Tisco) take their 
Keplerian values from Eq. (C.2) evaluated at the ISCO radius, 


while the radial component is fixed by unit-normalization: 
2 2 
2 ISCO ISCO 
u, = gr-[-1-2"( ) — aM (ug ) |. 
This defines a circular equatorial flow at every radius outside the 


event horizon ry: its four-velocity is given by Eqs. (C.1)-(C.2) 
for r > sco, and by Eq. (C.3) for ry < r < rsco. 


(C.4) 


Appendix C.1.2: Keplerian thick disk 


We now wish to find a natural way of extending the equatorial 
Keplerian disk to a geometrically thick configuration. 

For cylindrical radii outside the ISCO (op > 6), we simply 
define the specific angular momentum f to take the same value 
at every height z above the equator as it does in the equatorial 
Keplerian disk model. 
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That is, for 9 > 6, we impose the axially symmetric profile 


pr? 
&(p,2) = —>, (C.5) 
p-2 
resulting in a circular four-velocity of the type (C.1), 
f pr! 
gence dx" = -u,(p, z) (- dt + ue). (C.6) 
p-2 
whose unit-normalization fixes the time component to be 
p -1/2 
= tt 
—u,(P, Z) = (-« a (p 2 md ’ (C.7) 


where we used the fact that g# = (rsin@)* = p~?. Note that ¢” 
depends on z via the spherical radius r = +/p? + z2. It is easy to 
check that this quantity is well-defined for p > 6 (i.e., that the 
expression in parentheses remains positive). 

For cylindrical radii inside the ISCO (op < 6), a natural ex- 
tension of the Cunningham (1975) prescription is to require that, 
at any given height z, the components u, and ug keep their ISCO 
values for all radii p < 6. That is, 


ythick:<ISCO qyH — ulSCO(z) dt + u,(p,z) dr + up °(2) dd, 


ft (C.8) 


where the (t, ¢) components us, = urg(P = pisco,Z) are 
obtained by evaluating Eqs. (C.6)-(C.7) at the ISCO radius 
pisco = 6, while the radial component is again fixed by unit- 
normalization: 


[1,(0, 22 = fer (-1 — gi [8] — gi [us%]’). (C9) 


A simple numerical investigation of this relation reveals that it 
quickly turns negative as the height z above the disk increases, so 
that this four-velocity is not well-defined. We must therefore find 
a new prescription that produces a well-defined flow everywhere. 


Appendix C.1.3: General thick disk 


We will start from a general circular four-velocity of type (C.1), 


us (p, z) dx" = —u,(p, z) (— dt + €(p) dg), 


which assumes that the specific angular momentum depends 

only on cylindrical radius. The function ¢(p) is undetermined at 

this stage, but it is constrained by the requirement that this four- 

velocity remain well-defined everywhere outside the horizon. 
Its angular velocity is 


(C.10) 


g by £ 2 
G2. es = ;(1- ) (C.11) 
u stu r 
and unit-normalization fixes its time component to be 
- 1-2/r 
-u, = (ett — of) 1? — C.12 
u, = (-g" - 20?) eee (C.12) 


which is well-defined outside the horizon (r > 2) provided that 


Q€ <1, (C.13) 
or equivalently [by Eq. (C.11)], 

f 2 

S(i- *) « 1 (C.14) 
p 


Since p = rsin@, we have p — 2sin@ > 0 outside the horizon, 
where the condition (C.14) becomes 


pl? 


Vp - 2sind 


We have thus far imposed no conditions on the function €(p). 
At this point, we will assume that it converges to the Newtonian 
profile €Newton = pl! 2 at spatial infinity. We also want to consider 
a family of profiles that (i) includes the Keplerian profile (C.2) 
as a special case and (11) such that the constraint (C.15) leads to a 
simple condition. We are thus naturally led to consider the ansatz 


|&(p)| < (C.15) 


pr? 


pta 


la(p) = (C.16) 


where @ is areal constant such that @ = 2 recovers the Keplerian 
profile (C.2). However, we want our profile to be defined at every 
spacetime point and to remain free of singularities. We will thus 
require that a > 0, so that the denominator may never vanish. 
This also ensures that € > 0, and reduces condition (C.15) to 

(op +a) —p+2sind > 0. (C.17) 
Demanding that the discriminant of this quadratic equation be 
strictly negative results us in a stronger condition 


a>-, 


ri (C.18) 


which is the final condition for our angular momentum profile 
to be well-defined everywhere outside the horizon. Following 
Gold et al. (2020), we consider a = | in this article. We have 
numerically checked that this choice leads to a well-defined four- 
velocity even when the spin is nonzero. 

Finally, we note that there is still a locus of spacetime where 
our chosen profile (C.16) produces singular behavior: the axis 
p = 0. Taking the limit p — 0 keeping z > 2 fixed results in 


1/2 1/2 
v=(1-2} ; pti 2) ; 
r 


7 
so the ¢ component of the four-velocity, as well as the angular 
velocity Q, diverge in the limit p — 0. Indeed, Eq. (C.11) shows 
that the angular velocity grows like Q ~ fp~? ~ p7'/? near the 
axis for the Keplerian £, which behaves as € ~ p*/*. Hence, any 
family €(p) that includes the Keplerian profile as a member must 
have a divergent © near p = 0. However, since the orbiting mod- 
els that we consider in this paper always have vanishingly small 
emission close to the axis, this pathology never causes problems. 


(C.19) 


Appendix C.2: Infalling geodesic motion 


In this section, we work in the Kerr metric and consider 


Und 0, =u' 0, +u' 0, + u? Og, (C.20a) 


u,=—1, Uy = 0. (C.20b) 
This four-velocity describes particles that fall into the black hole 
from spatial infinity, where they start out with vanishing velocity 
and angular momentum. The (f, ¢) contravariant components are 


tt 


ui = guy =-§g, u? = gu, = —"*. (C.21) 
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Note that u* must be nonzero when the black hole is rotating 
(frame-dragging prevents purely radial geodesic infall). Unit- 
normalization fixes the r covariant component to be 


=| = tt 
uy = — eee 


= (C.22) 
& 


where we picked the negative sign of the square root to describe 
infall, and then the r contravariant component is 


u" = gu, = — y(-1 — 8") 9”. 


These formulas are easier to interpret in the Schwarzschild 
case, where they simplify to 


ae 2 
v=(1-2) : ul =—4l/-, 


r 


(C.23) 


(C.24) 


This results in a purely radial velocity (since ug = 0 at zero spin) 


dr ul 
e, = 


— Fr r= — ro 2 
veve=s “ie (C.25) 


where e, the unit spacelike vector e, = 0,/ /g;, and 


: V2 2) 2 
v=o = — y =, 
r r r 


which at spatial infinity reduces to the Newtonian formula for 
velocity derived from classical conservation of energy: 


(C.26) 


;WF-==0, (C.27) 


Appendix D: Precision of ray tracing computations 


In this section, we examine the strongly lensed null geodesics in 
Kerr, i.e., those that explore the shell of spherical photon orbits 
before escaping to produce high-order images. Their trajectories 
are very highly bent, as illustrated in Fig. D.1. 

Kerr null geodesics obey (e.g., Gralla & Lupsasca 2020a,b) 


me de & do 
[= ———— = ——. = Gs, D.1 
7 +, VR) A +,VOO ia 


where (r) and @(6) are radial and angular geodesic potentials, 
and the integrals are evaluated along the entire trajectory with the 
signs +, flipping at turning points of the radial/angular motion, 
which correspond to zeros of their respective potentials. 

We consider a geodesic traveling from a source at (r,, @,) to 
an observer at (ry, 6). To assess the precision of our geodesic 
integrator, we plot in Fig. D.2 the numerical evolution along the 
geodesic of the quantity |J, — Gel, as well as of the conserved 
energy E = —p, and azimuthal angular momentum L = pg. 

These quantities should remain constant along the geodesic, 
and the figure shows that their error is consistently below 107!, 
except for a peak of 107° in |J,, — Go| close to an angular turning 
point. We have also checked that the radii of the three equatorial 
crossings match the analytical formula derived by Gralla & Lup- 
sasca (2020b), which gives the radius of the equatorial crossing 
after m angular turns in terms of the Jacobi elliptic sine sn(x|k): 


2(1 (m) 
im _ Tarai Para st (3 Vrairateq — Fol k) 2) 


eq 2(1 (m) Flk 
r31 — fa, SN“ | 5 Y931142Teq — To 
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Fig. D.1. A null geodesic of order n = 2 (two @ turning points), in blue. 
The black disk is the event horizon and the yellow lines indicate the 
angular turning points 6_ and 6, = a — 6_, which are computed from 
the analytical expression for 6_ as a function of black hole spin and the 
geodesic constants of motion (see, e.g., Gralla & Lupsasca 2020a). The 
three equatorial crossings of the geodesic are labeled by the black dots. 
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Fig. D.2. Evolution of the difference between the initial and current 
values of the geodesic constants E = —p, (blue) and L = pg (orange), 
as well as the quantity |/, — Ge| (green), as a function of the integration 
step. The peak of |J — Go| happens close to an angular turning point 
where the geodesic evolves in the Kerr spherical null geodesics region. 


Here, rj; = rj — rj with {r), 72,73, r4} the four roots of the quartic 
potential R(r) (with exact expressions given in Eq. (A8) of Gralla 
& Lupsasca 2020a), and Ff, is the elliptic integral of the first kind 
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The Mino time elapsed up to the m' equatorial crossing is 


2mK — sign (P) F, 


(m) 

e — ? D.4 
Teg EE (D.4) 
where p® is the polar photon momentum at the observer, 

2 
7 > 70 _l ntra 
Us = Ag t [AG + 5, bo= 5(1- 2) (D.5) 


are the zeros u = cos” @ of @(u) in terms of the energy-rescaled 

angular momentum A = L/E and Carter constant 7 = Q/E7, and 
._ [Cos 

F,= F aresin( 2 


at), cer( 


u_ 2 

are an elliptic integral of the first kind and its completion K. 

We found that the three numerical and analytical values of 
agree to within 10-*—10-° M. Note that the equatorial cross- 
ing labeled 7°47 in Fig. D.1 occurs after the backward ray traced 
geodesic has encountered two angular turning points; this shows 
that a very high level of numerical precision is maintained de- 
spite the spike in the error |J, — Go| at that crossing. Finally, we 
have compared the value of the radial turning point rq (i.e., the 
minimum value of the radial trajectory) to the analytical expres- 
sion given in Eq. (A8d) of Gralla & Lupsasca (2020a), and found 
an error < 10~°M. All these results lead us to conclude that our 
null geodesic integrator is highly accurate, even within the pho- 
ton shell of bound spherical photon orbits. 

Besides the geodesic integration, Gyoro also integrates 
the radiative transfer equation by discretizing the part of the 
geodesic that intersects the disk, using a constant step of size 
6 = 0.1M. We have checked that reducing this step size by a 
factor of 10 changes the observed specific intensity by less than 
0.1%, which is sufficient given that our synchrotron emissivity 
prescription is only at the percent level of precision (see App. A). 
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Appendix E: Image orders 


Null geodesics in Kerr can be conveniently parametrized by the 
so-called Mino time (Gralla & Lupsasca 2020a,b) 


t=1,=Go, (E.1) 


with J and Gg are defined in Eq. (D.1). Before the first angular 
turning point (i.e., between the points S and 7, in Fig. 2), the 
Mino time elapsed along the geodesic as it travels backwards 
from the observer at 6, to 0 = @; is (Gralla & Lupsasca 2020b) 


7(6) = sign (p$) (Gp - Gi). (E.2) 
where p? is the polar photon momentum at 6 = @,, while 
1 7 
Go=- F (aresin( “*), (E.3) 
av—-u_ uy }| u- 


and G, and Gj, are evaluated at the Gp evaluated at 6 = 6, and 
0 = Os, respectively. We can determine the Mino time tT, elapsed 
at the first angular turning point by setting 


6, = arccos (+ vitz) : 


while the Mino time elapsed between two successive angular is 


+ = sign (P°) ; (E.4) 


(E.5) 


The order of any null geodesic is then defined as 


e n= 0: between the observer and tT = 71, 
e n= 1: between Tt = T; and t = 7] + Art, 
e n= 2: between tT = T; + At and t = 7] + 2AtT. 


The integration is cut off past n = 2, so as not to pollute the 
image by unresolved n > 2 features. 

We are thus able to provide not only the full image (con- 
taining all layers n < 2), but also individual images of each layer 
n € {0, 1,2} consisting of photons loaded onto each null geodesic 
at the corresponding order. 
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